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i.  Objectives: 

This  report  investigates  the  use  of  Bi-spectral  related 
techniques  to  extract  detection/classification  clues  from  time- 
frequency  representations. 

Earlier  results  have  indicated  that  1.5-D  spectral  techniques  (a 
degenerate  version  of  the  Bi-spectrum)  has  potential  Signal  to 
Noise  Ratio  (SNR)  gain  over  conventional  techniques.  This  is 
partially  due  to  the  rejection  of  Gaussian  like  perturbations  by 
the  cumulant  based  techniques.  The  third  order  moment  for  Gaussian 
zero  mean  random  processes  is  essentially  zero.  It  is  assumed  that 
i)  the  independence  of  IPS  from  analytic  signal  representations 
and  ii)  the  Gaussian  character  of  the  noise  causes  minimal 
distortion  of  the  output  variables,  hence  permit  signal 
characterization  at  low  to  moderate  SNR's. 

The  performance  is  demonstrated  on  i)  synthetic  signals  and  ii)  on 
segment (s)  of  the  NOSC  data  set. 


1.  Introduction: 

This  report  examines  the  instantaneous  power  spectrum  (IPS) 
and  the  cumulant  based  version  of  IPS  (CUM)  to  address  the  question 
whether  or  not  these  techniques  offer  possible  improvements 
relative  to  conventional  spectral  processing  (i.e.  spectrogram, 
Lofargram,  periodogram) . 

Section  2  provides  the  mathematical  definition  for  IPS,  CUM, 
and  LOFAR  while  section  3  allows  a  theoretical  comparison  of  these 
techniques.  The  comparison  is  performed  using  a  deflection 
criterion  (i.e.  maximum  signal  magnitude  versus  standard  deviation 
of  the  detector  output  under  noise  only  conditions) . 

Section  4  contains  simulation  results,  while  section  5 
provides  actual  SONAR  data  results. 

Section  6,  and  7  contain  the  conclusion  and  program  listings, 
respectively.  Section  8  provides  the  list  of  references. 
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2.  IPS,  CUM  and  LOFAR 


2.1  IPS: 

IPS  is  a  member  of  the  Cohen's  class  [1],  and  is  also  known  as 
the  real  part  of  the  Rihaczec  distribution.  Our  particular 
implementation  differs  in  that  filtering  is  employed  as  the  data  is 
processed  (equation  1) .  Earlier  results  have  shown  that  in  contrast 
to  the  Wigner-Ville  Distribution  (WVD)  no  spectral  cross  modulation 
terms  are  generated.  However,  it  has  wider  spectral  peaks  and 
superimposed  spectral  auto  modulation  terms  [2].  These  auto  terms 
do  not  interfere  with  the  interpretation  of  T-F  surfaces,  but  can 
help  in  the  detection  of  a  weak  stable  spectral  line  component.  The 
digital  implementation  of  IPS  is  given  by 

IPS{n,k)  {x(n)x‘ in-m)  *  x'  (n)  x(n+m)  )  w(m)  exp-j2Ttkm/N 

eq.  1 

where  k  =  frequency  index,  n  =  time  index,  N  =  the  length  of  data 
used,  m  =  shift  parameter  and  w(  )  =  lag  window. 

Some  earlier  results  are  published  in  [2,3,4]. 


2.2  CUM: 

In  general  the  Bi-spectrum  is  defined  as 


0)1,02)  =  ^ 


-j 


where  the  cumulant  is  given  by  [5] 


=  E{X‘  {n)X{ml^)X[n^l^))  . 
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A  degenerate  version  of  the  cumulant  is  given  by 


Cj^(0,22)  =  £(XUn)X(i2)X(n^22)  )  =  E(lX(n)  X{n*22)  )  . 


When  replacing  the  expected  value  with  the  instantaneous  value 
the  degenerate  Bi-spectrum,  also  called  the  Ih  D  spectrum,  becomes 

\^x(n+2) 

The  cumulant  version  of  IPS  is  given  by 
CUM(n,k)  =  —  (\x{n)  \^x*  {n-m)  +  |x(n)  Px-(n+7n)  )  w(m)  ex.p~j2nkm/N 

eq.  2 

Cross  terms  at  locations  where  the  true  spectrum  has  no  support 
are  not  created.  Also  contrary  to  the  conventional  Bi-spectrum  none 
of  the  spectrally  related  components  are  suppressed. 

In  part  of  the  work  in  [6],  the  approach  of  Petropulu  [7]  which 
uses  an  instantaneous  higher  order  moment  slice  was  utilized.  This 
approach  works  well,  but  does  not  permit  separation  of  spectral 
components  having  the  same  spectral  dynamics.  That  is,  all 
stationary  spectral  components  (i.e.  stable  line  components)  map  to 
the  same  detection  cell,  making  this  technique  useless  in 
scenarios  with  several  similarly  behaving  spectral  components. 


2.3  LOFAR: 

The  LOFAR  algorithm  uses  a  Hamming  windowed  Fast  Fourier 
Transform  (FFT)  and  displays  the  magnitude  of  the  transform. 


5 


To  keep  T-F  dimensions  relative  small,  no  zero  padding  is  used. 


P^(n,k)  =  I  x{w)  | 


eq.  3 

The  IPS,  CUM,  and  LOFAR  surfaces  are  averaged  (in  magnitude)  over 
some  segments  of  the  time  axis  leaving  only  a  display  of  spectral 
magnitude  versus  frequency.  For  each  segment  one  spectral  line  is 
created  resulting  again  in  a  time  frequency  display.  The  averaging 
technique  is  also  known  as  power  or  incoherent  averaging. 

3.  Theoretical  results  for  averaged  surfaces. 

The  analysis  in  this  section  is  done  by  disregarding  any 
window  effects  and  assuming  that  under  the  noise  only  hypothesis 
(Hq)  Gaussian  white  noise  and  under  the  signal  only  hypothesis  (H,) 
a  deterministic  signal  of  the  foirm  x(n)=A  cos (27rkQn/N)  is  present. 
The  value  of  kg,  n,  and  N  are  assumed  to  be  integers  which  ensures 
that  the  spectrum  peaks  at  a  bin  center  of  an  N-point  Fourier 
transform.  In  all  detection  schemes  (averaged  IPS,  averaged  CUM, 
and  averaged  LOFAR)  magnitudes  are  averaged  since  they  are 
numerically  more  conveniently  computed  then  the  magnitude  squared 
values. 

The  comparison  is  done  for  normalized  versions  of 
I(k)  =  llPS(n,k) I  , 

C(k)  =  |cUM(n,k)l  ,  and 

P(k)  =  [windowed  Fourier  transform  {x(n,,k))|  ;  where  n,  is  the 

time  counted  from  reference  point  n. 

In  what  follows  the  time  (n)  and  frequency  (k)  dependency  is 
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usually  suppressed. 


3.1  Averaged  IPS. 


IPS(n,k) 

real 

1 - 

- >  '  ' 

.  1  1 

1 

VM 

y  z  w 

where  x(n)  ~  N(0,a^)  i.i.d.  and  M  is  the  appropriate  number  of 
degrees  of  freedom. 


Under  Kg  : 

Under  noise  only  condition  (Hg)  the  random  variable  Y  denoted  by  Yg 
has  the  following  moments  [6]: 

E{Yg)=  a2  ,  E{Yg)2  =  [(N+5)/2]  ,  var(Yg)  =  [(N+3)/2]  a"  . 

For  tne  transformation  Z  =  SQRT(Y^)  , 

the  resulting  probability  density  function  is  given  by 
f2(z)  =  (fY(z)  +  fy(-z))  U(z)  . 

Making  the  assumption  that  for  large  enough  number  of  data  points 
the  random  variable  Y  is  approximately  Gaussian  distributed  we  can 
derive  the  moments  for  Z  as 
E(Zg)  =  SQRT{  (N+3)/167r} 

E(Zg2)  =  E{Yg)2  =  ((N+5)/2) 

var(Zg)  =  t(N+5)/2  -({(N+3)  }/167r)  ] 

var(Wg)  «  i/M  var(Zg)  (for  large  N) 
a  (1/M) 0.48  N  . 

A  typical  surface  has  about  N/3  degrees  of  freedom  [6],  providing 
a  standard  deviation  (under  the  Hg  hypothesis)  of 
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a(Wo)  »  1.2  . 

Under  (at  spectraj.  location  ±  k^) : 

Z,  =  IyJ  =  (NA^)  cos^CZjrkjjn/N) 

W,  =  (Na2)/M  cos2(2»rkon/N) 
a  (Na2)/2. 

The  deflection  Dp^,  for  IPS  becomes 

Dpj,  =  W,/a(Wo)  =  (NAV2)/(1.2  a^)  =  (N/2 . 4 )  (A^a^)  - 

3.2  Averaged  CUM. 


x(n) - > 

CUM 

- 5 

real 

1  1 

1  1 

- - J//M  E„ 

c 

Y 

z 

Under  Hg: 

E(Cg)  =  0  +  jo 
E(Cg)2  =  var(Co) 

=  (6N+54)/4  (a^)^  =  var(Yo)  ;  [6] 
(since  real  part  of  cum  is  zero) . 

Zq  =  I  ^0  I 

f^{z)=  2/(27ra/)’'2  exp- (zV  (2a/)  )  U(z) 
E(Zg)  =  SQRT(2/7r) 

E(Zo)2  =  3/2  a/ 

var(Zo)  =  (3/2  -  2/7r)  a/ 

E(Wo)  =  E(Zo)  =  SQRT(2/7r)  a^ 

E(Wg)2  =  1/M  {E(Zo)2  +(M-1)e2(Zo)  } 
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var(WQ)  =  1/M  {  (3/2-2/7r)  (6N+54)/4) 
aCWjj)  =  SQRT{1/M  (3/2-2/?r)  (6N+54)  }/2  . 

A  typical  surface  has  about  N/3  degrees  of  freedom  [6]  ,  so  the 

standard  deviation  becomes 

<7(Wo)  «  SQRT{3/N  (3/2-2/7r)  (6N+54)  }/2  . 

Under  H,  (at  spectral  location  ±  )  : 

x(n)  =  A  cos  (27rkQn/N) 

Y,  =  NaV2  cos^(27rkn/N)  @  k  =  ±  k^ 

Z,  =  NaV2  !  cos^(27rkn/N)  1  §  k  =  ±  k^ 

From  numerical  evaluations  we  know  that  the  summation  is  lower 
bounded  by 

—  (cosM27tln/M)  )  I  >  0.424  . 

Hence  Z,  >  0.424  NaV2. 

The  deflection  becomes 

=  (0.424  NAV2  )SQRT(M)/{SQRT[  (3/2-2/7r)  (6N+54)/4]  aM 

~  0.1075  N  (A/a)^  . 


3.3  LOFAR: 


x(n) - 


X(k) 


Under  H,,: 

x(k)=  x^(k)  +  j  x,(k) 
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r**N-1 


x(in)  exp-j  (2jrkin/N)  ;  where  n  is  the  reference  time 


E(X,(k))=  E(X^)=  E(X^(k))=  E(X.)=  0 
E(X^)2  =  var(X^)  =  var(X^)  =  N/2 

Z  =  SQRT(X^2  +  Xj2  ) 

E(Zo)  =SQRT(Tr/2)  SQRT(N/2)  a 

var(Zo)  =  (2-7r/2)a2  N/2 

E(Wq)  =  E(Z(j)  =SQRT(jr/2)  SQRT(N/2)  a 

E(Wo)2  =i/M  {E(Zo)2  +  (M-1)e2(Zo)) 

var(Wo)  =  1/M  var(Zo)  =  (1/M)  (2-Jr/2)  (N/2) 

a(Wo)  =SQRT{(1/M)  (2-jr/2)  (N/2)}  a 

Under  (at  spectral  location  ±  k^  ) : 

X,  (k)  =A  N/2  §  k  =  ±  kg 

Z,  =  A  N/2  @  k  =  ±  k(j 

E(W^)  =  E(Zi)  *  A  N/2 
The  deflection  becomes 

=  (AN/2)  /{  (1/M)  (2-7r/2)a2(N/2) 

=  n’/^  [3  A^  /(  (4-7r)a2)  ]  . 

3.4  Deflection  ratio  comparison: 

We  can  now  plot  the  deflection  ratio,  that  is 
determine  which  detector  is  more  advantageous 
conditions. 

Comparing  D,pg  and  : 

^IPS  ^  ^LOFAR 

simplifies  to 
N  >  20.13  (oVa^) 

Comparing  and  : 
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compare  and 
under  which 


^CUM  ^  *^LOFAR 

simplifies  to 
N  >  302.4  (oVa^)  . 

We  notice  that  the  processing  length  variable  N,  for  IPS  relative 
to  LOFAR  has  a  (SNR)'^  and  for  CUM  relative  to  LOFAR  has  a  (SNR)'^ 
dependency. 


4.  Simulation  results  for  sinusoid  in  white  Gaussian  noise. 

Ten  test  sets,  of  length  4096,  consisting  of  sinusoidal 
signals  embedded  in  white  Gaussian  noise  are  processed  by  IPS,  CUM 
and  LOFAR.  Each  test  set  contains  seven  (7)  sinusoids  whose  SNR 
goes  from  -3  dB  to  -21  dB  in  3  dB  steps.  The  decrease  is  monotonic 
relative  to  the  spectral  locations  so  that  the  SNR  decreases  in  3 
dB  steps  as  the  frequency  increases.  The  sinusoidal  frequencies 
are  fixed  in  all  test  sets  at  digital  frequencies: 

{w  )  =  {11%,  17%,  29%,  37%,  59%,  77%,  111%)  *27r/256  .  The  phases 
within  each  set  and  over  all  set  are  independently  distributed 
(i.e.  no  harmonic  relationship  between  different  spectral 
components  and  different  realizations) ,  We  note  that  none  of  the 
sinusoids  has  an  integer  number  of  cycles  over  4096  data  points 
ensuring  all  signals  will  have  side  lobes  if  processed  by  an  FFT  of 
size  4096  or  less.  To  demonstrate  potential  gain  the  deflection 
ratio  as  defined  in  section  3  is  computed  for  each  sinusoid  and 
averaged  over  the  10  realizations. 

Figures  1,  2,  3,  4,  5  and  6  are  plots  of  the  outcomes  of  the 
simulation  for  128,  256,  512,  1024,  2048  and  4096  data  points. 
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respectively.  The  top  part  of  each  figure  superimposes  all  of  the 
10  realizations  while  the  bottom  part  shows  the  average  of  the  10 
realizations.  We  note  for  example,  that  the  weakest  sinusoid  (bin 
223-224  of  512  point  data  length)  clearly  shows  up  in  the  CUM  and 
IPS  representation  (figure  3. A)  but  is  not  detectable  in  the  LOFAR 
representation  (figure  3.B).  This  illustrates  that  the  alternative 
representations  (i.e.  averaged  TPS  or  averaged  CUM)  can  complement 
the  traditional  (averaged  LOFAR)  representation.  Figures  7-12 
provide  plots  of  the  experimental  deflection  ratios  versus  SNR  in 
the  top  half  of  the  figures.  The  solid,  dashed,  and  dotted  lines 
represent  the  experimental  deflection  for  CUM,  IPS,  and  LOFAR, 
respectively.  The  bottom  halves  are  plots  of  the  difference  of 
deflection  ratios  (Djpg  -  )  .  Figure  13  attempts  to  reconciliate 
theoretical  bounds  from  section  3  with  the  experimental  evidence 
from  this  section.  The  straight  solid  lines  demarcate  the 
theoretical  transitions  points  at  which  CUM  works  better  than  LOFAR 
(upper  straight  line)  and  where  IPS  works  better  than  LOFAR  (lower 
straight  line) .  We  mention  again  that  no  window  action  was 
accounted  for  in  the  derivation  of  these  lines.  Superimposed  are 
(zigzag  lines)  the  results  taken  from  figures  7  trough  12.  The 
upper  zigzag  line  shows  the  demarcation  line  where  the  experimental 
evidence  indicates  IPS  will  do  better  than  CUM  (in  terms  of 
deflection  ratios) .  The  lower  zigzag  line  consists  of  the  points, 
as  a  function  of  SNR  and  processing  length,  where  the  experimental 
performance  of  all  three  techniques  (IPS,  CUM  and  LOFAR)  tends  to 
coincide.  Above  this  line  the  CUM  and  IPS  technique  outperform  the 
LOFAR  approach  with  improvements  increasing  as  the  input  SNR  increases. 
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Fig.  l.A 
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5.  Processing  Results  using  N08C  Data  Base. 

In  this  section  the  data  set  supplied  by  NOSC  is  processed. 
Figure  14  is  a  LOFARGRAM  of  the  data  set  r2A_14  which  is  4,128,769 
data  points  long.  The  time  origin  for  all  grains  is  chosen  to  be  at 
the  top  of  each  figure.  Figure  15. A  and  15. B  are  LOFARGRAMS 
produced  at  NFS  to  verify  the  fidelity  of  the  data  transfer  from 
NOSC  to  NFS.  To  see  representative  spectral  plots  we  selected  the 
most  interesting  segment,  that  is  from  time  sample  1,024,001  to 
1,228,800.  This  corresponds  to  the  region  508  to  408  along  the  time 
(vertical)  axis  of  figure  15. A.  Initial  processing  using  IFS  and 
especially  CUM  showed  heavy  modulation  of  the  time  frequency 
displays.  This  is  caused  by  the  low  frequency  spikes  observable  on 
figure  14  and  15. A  (i.e.  strong  DC  like  pulsations  (about  5800 
sample  points  apart) .  In  the  top  part  of  figure  16  which  displays 
the  time  series,  used  in  the  experiments,  these  modulation  spikes 
are  clearly  seen.  To  minimize  the  unwanted  amplitude  modulation 
effects  on  the  time  frequency  plots,  the  data  is  soft  limited  by 
scaling  all  values  larger  than  +90  by  dividing  them  by  4  and  by 
scaling  all  values  less  than  -90  by  dividing  them  by  2.  The 
resulting  time  series  is  plotted  in  the  bottom  part  of  figure  16. 

The  contour  plots  of  the  remainder  of  this  section  show  now 
little  amplitude  modulation  allowing  the  approximation  to  gray 
scale  via  contour  plots  to  work.  On  a  color  monitor,  with  a 
sufficient  number  of  color  levels  of  quantization,  the  modulation 
may  not  be  bothersome.  Figure  17  displays  49  traces  (along  the  y- 
axis)  which  corresponds  to  49  individual  estimates  of  the  spectral 
density.  Closer  examination  reveals  some  details  on  the  IFS  plot 
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not  so  easily  seen  on  the  lower  LOFARGRAM  plot.  In  particular  we 
refer  to  3  areas: 

1)  tick  mark  38-48  (vertical  axis)  about  bin  700  (horizontal 

axis)  , 

2)  tick  mark  1-33  (vertical  axis)  about  bin  850  (horizontal 

axis) ,  and 

3)  tick  mark  1-20  (vertical  axis)  about  bin  1300  (horizontal 
axis)  . 

A  word  of  caution  however  when  using  the  simple  minded  contour 
(MATLAB)  plots:  peaks  of  identical  height  but  of  different  width 
will  appear  differently.  In  appendix  A  the  effects  of  choosing 
different  contour  levels  is  illustrated.  However  for  all  plots  the 
best  choice  of  level  was  manually  selected  by  examining  several 
plots  differing  in  the  number  of  levels  selected  and  then  retaining 
the  ono  with  the  best  features. 

Figure  18. A  and  18. B  are  displays  of  IPS,  CUM  and  LOFAR  with 
the  lines  created  512  time  data  points  apart.  Figure  18. A  (as  well 
as  19. A  and  20. A)  displays  spectral  location  51  through  1500  while 
figure  18. B  (as  well  as  19. B  and  20. B)  displays  spectral  location 
351  through  1500.  This  allows  better  interpretation  of  the 
different  spectral  regions. 

Figures  18,  19  and  20  have  the  same  processing  parameters, 

they  differ  in  the  number  of  spectral  lines  computed.  Figure  18 
uses  a  shift  of  512  (i.e.  8:1  overlap)  resulting  in  390  lines, 

figure  19  uses  a  shift  of  1024  (i.e.  4:1  overlap)  resulting  in  197 
lines,  while  figure  20  uses  a  shift  of  4096  (i.e.  no  overlap) 

resulting  in  49  spectral  lines. 
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Figure  21  illustrates  performance  as  the  processing  length  is 
increased  to  8102  (shift  8192,  no  overlap)  providing  25  spectral 
lines  (bin  51  through  3300  are  displayed) . 

6.  Conclusion  and  recommendation. 

As  illustrated  in  section  3  and  4  processing  gain  can  be 
obtained  for  IPS  and/or  CUM  relative  to  the  LOFAR  processing 
technique.  This  is  a  function  of  signal  and  processing  parameters. 
It  may  allow  glimpses  of  targets  which  otherwise  maybe  hidden. 
However,  the  processing  gain  when  achievable  is  bought  at  the 
increase  in  processing  cost.  IPS  and  CUM  are  actually  tools 
designed  more  for  transient  type  signals  but  as  shown  can  be  used 
to  obtain  gain  on  stable  narrow  band  signals.  Simplified 
processing  was  not  attempted  but  deserve  future  examination.  The 
analysis  in  section  3  is  done  totally  disregarding  lag  and  data 
windows  which  typically  are  used.  Future  analysis  should  address 
the  influence  of  the  window  function  on  processing.  Extensions  to 
extract  phase  information  is  not  attempted  but  should  be  pursued. 
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7.  Program  Listings: 

All  algorithms  are  implemented  (including  their  hard  copy 
functions)  using  MATLAB  software  (The  Mathworks  Inc.) .  The  programs 
can  be  executed  on  MS-DOS  personal  computers  (i.e.  386  +  co¬ 
processor  or  486  based  IBM  compatible  PC's)  or  on  SUN  work 
stations.  Due  to  the  long  data  sequence  and  large  T-F  surfaces  most 
of  the  work  is  done  on  SUN  work  stations.  The  hard  copies  provided 
are  obtained  by  using  the  contour  plotting  function  of  MATLAB.  One 
note  of  caution  is  appropriate  when  interpreting  the  T-F  plots.  Due 
to  the  way  contour  lines  are  drawn,  equal  height  spectral  peaks  can 
appear  different  if  they  differ  in  T-F  spread,  that  is  a  narrow 
peak  will  be  darker  than  a  wide  peak  even  though  these  peaks  are  of 
the  same  height.  This  is  due  to  the  density  of  contour  lines  which 
is  higher  for  narrow  peaks  than  it  is  for  wide  peaks.  Better 
results  are  obtained  when  displaying  the  data  on  gray  scale  or 
color  scale  output  devices.  The  new  4.0  MATLAB  version  supports  the 
function  'image'  which  on  a  gray  scale  or  color  monitor  will 
improve  the  readability  of  LOFARGRAMs. 

The  color  output  (using  MATLAB  3.5i)  on  a  SUN  SPARC-2  work 
station  shows  superior  detail  when  compared  to  the  hard  copies 
supplied  in  this  report.  The  hard  copies  are  obtained  using  a 
screen  dump  routine  which  preserves  less  detail  than  a  postscript 
derived  output.  The  postscript  output  is  not  used  due  to  the  length 
of  time  involved  generating  it  on  the  processing  system. 
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>>  atypical  setup  for  obtaining  IPS  and  CUM  like  surfaces 
>>  ^for  i=start:finish 

>>  ^pCi ,:)=ipsavCrfC(i-1)*shft+1 : Ci-1)*shft+N),1 ,N,step); 

>>  ^qCi  ^:)=curaavCrfCO'-1)*shft+1 :  0'-1)*shft+N),1  ,N,step)j 

>>  %  N=  proc  length,  step=stepsize,  N/step=nuiiiber  of  terms  averaged  to 

>>  ^create  one  spectral  line,  shft  =  amount  of  data  shifted  to  create 

>>  ^the  next  surface  which  will  be  averaged  to  create  one  spectral 

>>  ^line  at  location  i 

>>  ^end 

>> 

>>  atypical  setup  for  obtaining  LOFAR  like  surfaces 
>>  ^for  i=start:finish 

>>  ^poCi ,:)=lofarinCrfCCi-1)*shft+1 : Ci-1)*shft+N),id,st,tl); 

>>  ^id=number  of  data  points  in  transform,  st=stepsize, 

>>  ^tl  =transform  length 

>>  %  Note:  rf  =  input  sequence, 

>>  %  Note:  in  all  run  executed  id=st=tl 
>>  %  So  the  overlap  was  controlled  with  shft 
>> 

>> 

XN  ■ 
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%function  [P, freqindex] =ipsav (data, wintype, winlen,  step)  ; 

%This  function  will  calculate  a  single  averaged  Instantan.  Power  Spectral  (IPS) 
%  line. 

%The  IPS  surface  characteristics  are  determined  by  the  selection  of  window 
%type  (wintype) ,  window  length  (winlen)  and  the  distance  that  the  window 
%is  moved  through  the  data  sequence  (step) . 

% 

%The  P  matrix  plots  only  positive  frequencies  (in  magnitude) .  The 
%outputs  timeindex  and  freqindex  can  be  used  in  plots  to  interpret  the 
%results . 

%The  inputs  are: 

%data  -  The  input  data  string  row  vector 
Iwintype:  '0'  Rectangular  Window 
%  '1'  Hcimming  Window 

%winlen  -  The  desired  width  the  window,  normally  half  of  the  siglen 
%step  -  Time  step  desired,  normally  '1'  or  a  multiple  of  '2' 

%See  also  IPSSURF,  IPSLOFAR 

function  [P, freqindex] =ipsav (data, wintype, winlen,  step) 

[datarows,  datacolumns] =si2e (data)  ; 
if  datarows  ~=1 
data=data . ' ; 

end 

siglen=length (data) ; 
if  wintype==0 

win=ones (winlen-1,  1)  ; 
elseif  wintype==l 
win=haraming (winlen-1)  ; 

end 

W= [win (winlen/2 :-l:l)]; 

>;=[  zeros  (1,  winlen)  data  zeros  (1,  winlen)  ] ; 

P=zeros (1, winlen/2) ; 

for  n=winlen+l : step : siglen +winlen-step+l 

Xm= [con j (x (n : -1 : n- (winlen/ 2- 1 ) ) ) . '  x (n :n+ (winlen/2- 1 ) ) . ' ] ; 

Xn=  [x  (n) ;  con:  (x  (n)  )  ] ; 
product= ( (Xm*Xn) . *W) . ' ; 

product= [product  0  con] (product (winlen/2 :-l :2) )]  ; 
ptemp=f f tshift (real ( , 5*f ft (product) ) )  ; 

P=P+abs  (ptemp  (winlen/2  +  1 :  wimen)  )  ; 
end 

[prow,  pcolumn] =size (P) ; 
freqindex= [ 0 :pcolumn-l ] ; 
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%function  [P,  freqindex, timeindex] =cuinav (data, wintype,  winlen,  step)  ; 

%This  function  will  calculate  the  1.5  D  Spectral  surface. 

%this  surface  consists  of  averaged  (shorter)  cumulant  based  surfaces 
%in  a  power  or  magnitude  averaged  sense 

%The  1.5  D  surface  characteristics  are  determined  by  the  selection  of  window 
%type  (wintype),  window  length  (winlen)  and  the  distance  that  the  window 
%is  moved  through  the  data  sequence  (step) . 

%The  surface  is  placed  sequentially  in  the  P  matrix  for  display. 

%The  P  matrix  plots  only  the  positive  half  of  the  spectral  plane.  The 
loutputs  timeindex  and  freqindex  can  be  used  in  plots  to  interpret  the 
%results . 

%The  inputs  are: 

%data  -  The  input  data  string  must  be  in  row  vector  form 
Iwintype:  '0'  Rectangular  Window 

%  '1'  Hamming  Window 

%winlen  -  The  desired  width  of  the  window,  normally  half  of  the  siglen 
%step  -  Time  step  desired,  normally  '1'  or  a  multiple  of  '2' 

%See  also  ONESURF,  ONELOFAR 


function  [P,  freqindex] =cumav (data, wintype,  winlen,  step) 

[datarows, datacolumns ] =size (data)  ; 

if  datarows~=l 

data=data . ' ; 

end 

siglen=length (data) ; 
if  wintype==0 

win=ones (winlen-1, 1) ; 
elseif  wintype==l 
win=hamming (winlen-1 ) ; 

end 

W= [win ( winlen /2 :-!:!)]; 

x= [zeros (1, winlen)  data  zeros  (1, winlen) ] ; 

P=zeros (1, winlen/2) ; 

for  n=winlen+l : step : siglen+winlen-step+1 
Xn=  [abs  (x  (n)  )  ^2  ;  abs  (x  (n)  )  ^2  ]  ; 

Xm= [con j (x (n : -1 : n- (winlen /2- 1 ) ) ) . '  x (n : n+ (winlen /2-1 )).']; 
product= ( (Xm*Xn) . *W) . '  ; 

product= [product  0  con^ (product (winlen/2 :-l :2) )]  ; 
ptemp=f ftshift (real ( . 5*f ft (product) ) ) ; 

P=P+abs (ptemp (winlen/2  +  1 : winlen) )  ; 

end 

[prow,pcolumn]=size  (P) ; 
f reqindex= [ 0 :pcolumn-l 1 ; 
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%  This  program  computes  the  LOFARGRAM  based  on  windowed  periodogramis  and  dispia%  ys  the  resu 
%id=no  of  points  in  xform,,  step  =  shift,  tlen=transfcrr.  length 
% [P] =lo farm (name, id, step, tlen) 
function  [P ] =lo farm (name, id, step, tlen) 
clear  pow 

%name=input ('  name  of  input  file  is  '  ) ; 

[drow,  dcolumn ] =size (name) ; 
if  drow'=l 

name = name . '  ; 

end 

len=length (name) ; 

%id=input('  No.  of  non-zero  data  points  in  transform  '); 

%step=input  ( '  shift  (step  size)  in  data  points  '); 

%tlen=input (' transform  length  '  )  ; 

i=fix ( (len-id) /step) +1; 
win=hamming (id) 
for  ic=l:i 

ppow=abs ( f ft (name ( 1+ (ic-1 ) »step : id+ (ic-1 ) *  step) . *win,  t len) )  ; 
pow(ic, : ) =ppow (1 :tlen/2) ; 
end 
P=pow; 

clear  name  drow  dcolumn  len  id  step  tlen  i  ic  win  ppow 
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9.  Appendix  A 


This  appendix  provides  3  figures  which  illustrate  the  effects  of 
level  selection  when  using  a  simple  (no-gray  scale  ,  no  color  ) 
printer.  All  three  figures  use  a  shift  of  4096  (i.e.  no  overlap  of 
data  between  successive  spectral  lines).  Fig  A.l,  A. 2  and  A. 3  show 
IPS,  CUM,  and  LOFAR  respectively.  The  number  of  quantization 
levels  is  different  on  each  figure.  The  label  on  each  figure  is 
self  explanatory. 
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